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INTRODUCTION 


According  to  Finlayson  and  Scriven^  it  is  not  variational  notation  or 
even  the  concept  of  a  varied  path  which  is  the  key  criterion  of  a  true 
variational  'principle',  but  rather  the  existence  of  a  functional  which  when 
varied  and  set  to  zero,  generates  the  governing  equations  and  constraints  for 
a  given  class  of  problems.  In  this  sense,  certain  fundamental  principles  of 
mechanics  such  as  d’Alembert's  Principle  do  not  truly  qualify  as  variational 
principles.  That  is  to  say,  these  mechanical  principles  or  'laws'  cannot  be 
posed  as  central  problems  of  the  calculus  of  variations.  On  the  other  hand 
there  are  others,  such  as  Hamilton's  principle  which  do  qualify  as  true 
variational  principles.  Yet  it  is  d'Alembert's  Principle  which  forms  a  basis 
for  all  analytical  mechanics* 2  and  it  follows,  therefore,  that  the  vanishing  of 
the  first  variation  of  some  functional  is  not  a  necessary  condition  for  the 
scalar  formulation  of  any  mechanics  problem  -  however  elegant  or  convenient 
this  may  be. 

Whether  a  true  variational  principle  or  a  more  fundamental  variational 
statement  is  used  to  obtain  a  numerical  solution  to  a  dynamics  problem,  an 
important  argument  is  that  well  established  laws  such  as  d'Alembert  s 
Principle  or  true  principles  such  as  Hamilton's,  are  physically  based  and 
avoid  the  arbitrariness  inherent  in  general  weighted  residual  methods  and 
contrived  variational  principles.  Only  variational  principles  which  are  also 
maximum  or  minimum  principles  appear  to  offer  any  advantage  for  obtaining 

finlayson,  B.  A.  and  Scriven,  L.  E.,  "On  the  Search  for  Variational 
Principles,"  Int .  J.  Heat  Mass  Transfer,  Vol.  10,  1967,  p.  799-821. 

2Lanczos ,  C.,  The  Variational  Principles  of  Mechanics,  3rd  Edition,  University 
of  Toronto  Press,  1966,  pp.  70-72. 
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approximate  solutions  -  mainly  through  their  ability  to  provide  bounds  on  the 
variational  integral.  Even  then  the  system  treated  must  be  positive-definite 
and  the  upper  and  lower  bounds  are  often  too  far  apart  to  be  of  practical 
value.  In  brief,  there  seems  to  be  little  point  in  contriving  a  variational 
principle  in  preference  to  a  variational  law  of  mechanics  despite  the  more 
primitive  status  of  the  latter.  Indeed  the  many  solutions  to  initial  value 
dynamics  problems  achieved  by  C.  Bailey^  by  applying  the  Ritz  method  to 
Hamilton’s  ’law  of  varying  action’  demonstrate  the  usefulness  of  variational 
formulations  not  qualifying  as  ’principles’.  Thus  motivated,  the  work  herein 
explains  the  numerical  difficulties  encountered  in  attempting  to  generalize 
Bailey’s  formulations  according  to  the  method  of  finite  elements. 

FINITE  ELEMENTS  IN  TIME 

The  many  solutions  achieved  by  C.  Bailey  were  generated  by  the  Ritz 
method^  using  a  power  series  approximation  in  which  globally  defined  poly¬ 
nomials  are  the  basis  functions.  Ultimately  the  length  of  interval  over  which 
solutions  may  be  generated  as  well  as  the  detail  to  be  provided  in  any  subin¬ 
terval  will  be  limited  by  the  degree  of  polynomial  used  as  a  basis.  The  pit- 
falls  of  using  higher  powered  polynomials  are  well  documented^  and  partially 
account  for  the  use  of  locally  (piecewise)  defined  basis  functions  (finite 
elements)  to  solve  problems  in  many  branches  of  mathematical  physics.  The 


^Bailey,  C.  D. ,  "The  Method  of  Ritz  Applied  to  the  Equation  of  Hamilton," 
Computer  Methods  in  Applied  Mechanics  and  Engineering,  7,  1976,  pp.  235-247. 
^Kantorovich,  L.  V.  ,  and  Krylor,  V.  I.,  Approximate  Methods  of  Higher 
Analysis,  Interscience  Publishers,  Inc.,  1964,  pp.  258-303. 

"’Conte,  S.  D.,  and  de  Boor,  C.,  Elementary  Numerical  Analysis:  An  Algorithmic 
Approach,  2nd  Edition,  McGraw  Hill,  1972,  pp.  231-233. 
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extraordinary  accuracy  and  simplicity  of  procedure  attained  by  Bailey, 
however,  is  not  to  be  understated. 

Apart  from  avoiding  the  problems  which  can  arise  when  higher  powered 
polynomials  are  employed  as  basis  functions,  finite  element  formulations  have 
other  advantages  when  used  to  solve  problems  in  continuum  mechanics.  Even 
though  the  principal  motivation  for  their  use  has  been  the  need  to  handle  com¬ 
plicated  boundary  shapes  (non-existent  in  the  time  domain)  finite  elements  are 
also  well  suited  to  handle  sudden  changes  in  load  functions,  extending  the  in¬ 
terval  of  solution  indefinitely  without  restart,  and  providing  great  detail  to 
the  solution  in  any  subinterval.  Thus  despite  the  reservations  expressed  by 
Zienkowicz ,6  the  extension  of  the  finite  element  method  to  the  solution  of 
transient  field  problems  is  well  motivated  and  was  first  reported  by  Argyris 
and  Sharpf^  and  later  by  Fried. ® 

Both  of  these  works  attempt  to  use  Hamilton’s  principle  as  a  starting 
point  for  the  finite  element  formulation  of  initial  value  problems.  As  will 
be  pointed  out  in  the  following  section,  this  cannot  be  accomplished  without 
some  logical  inconsistency  when  bringing  the  initial  data  into  the  formula¬ 
tion.  In  the  sequel  it  will  be  shown  that  the  use  of  Hamilton’s  ’law’,  rather 
than  Hamilton’s  ’principle’,  makes  possible  the  logical  incorporation  of  the 
initial  conditions  into  the  variational  formulation. 

^Zienkiewicz ,  0.  C. ,  The  Finite  Element  Method,  3rd  Edition,  McGraw-Hill, 

1977,  pp.  569-70. 

^Argyris,  J.  H. ,  and  Scharpf,  D.  W. ,  "Finite  Elements  in  Time  and  Space," 
Nuclear  Engineering  and  Design,  10,  1969,  456-464. 

®Fried,  I.,  "Finite-Element  Analysis  of  Time-Dependent  Phenomena,"  AIAA 
Journal,  7,  No.  6,  pp.  1170-1172. 
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HAMILTON’S  PRINCIPLE  -  A  CONSTRAINED  VARIATIONAL  PRINCIPLE 


The  following  equation  is  known  as  the  generalized  principle  of 
d ’ Alembert 

N 

l  ( F  i — P  i )  *  <5ri  =  0  ;  (')  =  3/3t  (1) 

i=l 

This  equation  applies  to  any  system  of  N-particles,  the  ith  particle  having  a 

position  r^,  a  momentum  P^  ,  and  subject  to  a  resultant  applied  force  F^. 

Under  the  assumption  that  the  virtual  work  of  the  applied  forces  is 

derivable  from  a  scalar  V,  a  time  integration  of  equation  (1)  leads  to 

Hamilton’s  law  of  varying  action: H 

1 2  N  1 2 

6/  (T-V)dt  -  l  mir^Sri]  =0  (2a) 

ti  l-l  ~  ~  tx 

T  is  the  kinetic  energy  of  the  system 


T  =  1/2  l 

i=l  "  ~ 

and  V  is  the  potential  energy  of  the  forces  impressed  on  the  N-particles.  The 
existence  of  V  makes  little  difference  as  far  as  numerical  calculations  are 
concerned.  In  the  event  V  does  not  exist,  equation  (2a)  can  be  written: 


1 2  -  N  •  t2 

/  (6T+S W)dt  -  l  m^  •  6ri]  =  0  (2b) 

tl  1-1  ~  ~  t i 


^Mierovitch,  L.,  Methods  of  Analytical  Dynamics,  McGraw-Hill,  1970,  p.  65. 
^Bailey,  C.  D. ,  "Application  of  Hamilton's  Law  of  Varying  Action,"  AIAA 
Journal,  Vol.  13,  No.  9,  pp.  1154-1157. 

^Hamilton,  W.  R. ,  "Second  Essay  on  a  General  Method  in  Dynamics," 

Philosophical  Transactions  of  the  Royal  Society  of  London,  1S35,  pp.  95-144. 
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The  bar  signifies  that  in  general  the  virtual  work  of  the  applied  forces 
cannot  be  derived  from  any  scalar  function  of  the  generalized  coordinates. 
Either  of  equations  (2)  can  be  used  as  a  basis  for  a  Ritz  approximation  to  a 
dynamics  problem. 

If  irj(ti)  and  6ri(t2)  vanish  in  equation  (2a),  the  result  is  Hamilton's 
principle : 


«/ 


c2 


(T-V)dt  =  0 


(3) 


tl 

Since  the  vanishing  of  the  displacement  variations  at  the  end  points  is 
not  the  only  means  by  which  the  partial  sum  in  equation  (2a)  may  vanish,  equa¬ 
tion  (3)  may  not  always  represent  Hamilton’s  principle  in  the  strict  sense. 
Should  equation  (3)  be  used  as  a  basis  for  the  numerical  solution  of  a  dynam¬ 
ics  problem  without  the  requirement  that  all  of  the  6r^  vanish  at  or  t2> 
zero  momentum  conditions  will  prevail  instead  as  natural  boundary  conditions 
on  those  displacements  whose  variations  are  free.  This  aspect  of  variational 
principles  is  covered  very  clearly  in  many  references  (cf.  ref.  12).  An 
observation  to  be  made  here  is  that  equation  (3)  corresponds  to  a  system  of 
boundary  va.lue  problems  -  not  initial  value  problems  -  since  the  partial  sum 
can  only  vanish  through  boundary  (end  point)  constraints  either  natural  or 
imposed.  Thus  equation  (3)  cannot,  with  complete  logic,  be  used  to  formulate 
any  system  of  initial  value  problems  of  dynamics.  The  introduction  of  initial 
data  has  in  fact  always  been  the  obstacle  preventing  the  use  of  Hamilton’s 


l^Couranc ,  R.  ,  "Variational  Methods  for  the  Solution  of  Problems  of 

Equilibrium  and  Vibrations,"  Bulletin  of  American  Mathematical  Society,  49, 
pp.  1-23. 
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principle  for  the  variational  formulation  of  initial  value  problems.  “ » ^ 

Since  equation  (3)  is  a  valid  physical  statement  of  mechanics  only  when 
the  boundary  constraints  are  such  that  the  partial  sum  vanishes,  it  is  proper 
to  refer  to  this  equation  as  a  'constrained  variational  principle'  as  opposed 
to  equations  (2)  which  are  unconstrained  variational  laws  of  mechanics,  suit¬ 
able  for  the  application  of  arbitrary  constraint  conditions. 


GLOBAL  AND  PIECEWISE  RITZ  APPROXIMATIONS 


Equations  (2)  and  (3)  differ  only  in  the  presence  or  absence  of  boundary 
terns.  For  the  case  of  a  single  particle  (N=l)  having  only  one  degree  of 
freedom  u(t),  the  Ritz  procedure  when  applied  to  either  of  equations  (2)  leads 
to  a  scalar  relation  of  the  form: 

6Ut[ (K-B)U-F]  =  0  (4) 

whereas  for  equation  (3) : 


6UT[ KU-F]  =  0  (5) 

Equations  (4)  and  (5)  are  assumed  to  derive  from  applying  the  Ritz  procedure 
whereby  the  displacement  function  u(t)  is  approximated  as: 

u(t)  =  a^(t)U  (6) 

The  relation  (6)  applies  to  the  entire  interval  of  solution  when  globally 
defined  basis  functions  are  used  or  to  a  particular  subinterval  thereof  when 
piecewise  functions  (finite  elements)  are  employed.  When  a  global  power  ser¬ 
ies  approximation  is  used  U  is  a  vector  of  generalized  coordinates,  the  first 


l^Tiersten,  H.  F.,  "Natural  Boundary  and  Initial  Conditions  From  a 

Modification  of  Hamilton's  Principle,"  J.  of  Math.  Physics,  Vol.  9,  No.  9, 
pp.  1445-1450. 

l^Gurtin,  M.  E. ,  "Variational  Principles  for  Linear  Elastodynamics , "  Archive 
Rati.  Mech.  Anal.  16,  34-50  (1964). 
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two  of  which  are  identifiable  as  u^)  and  u(tj).  The  'shape  function',  a(t), 
in  this  case  is  simply: 

aT(t)  =  [1  ,t  ,t2  .  ,tn]  ,  tj  <  t  <  t2  (7) 

If  piecewise  cubic  Hermite  polynomials  are  used  instead,  the  components  of  U 

are  local  values  of  u  and  u  defined  at  the  endpoints  of  a  particular  subinter¬ 
val  ,  and 

aT(t)  =  [2t3-3t2+1,  h( t 3  — 2t 2-f*T )  ,  3t2-2t3,  h(t3-x2)]  (8) 

where  t  =  t/h,  h  being  the  length  of  the  particular  subinterval.  Referring 
rirst  to  equation  (5),  it  is  noted  that  K  tends  to  be  singular  of  degeneracy 
one.  For  certain  simple  problems  K  may  compute  to  be  exactly  singular.  In 
general,  however,  K  will  only  become  singular  in  the  limit  as  the  number  of 
basis  functions  employed  in  the  Ritz  approximation  becomes  infinite.  The 
degeneracy  of  K  represents  the  possibility  that  neither  u(t!>  or  u(t2)  has 
been  specified.  That  is,  if  neither  6u(tl)  or  Su(t2)  vanishes,  then  mu  must 
vanish  at  both  endpoints  as  natural  boundary  conditions.  Under  these  condi¬ 
tions  u N t )  may  only  be  determined  to  within  an  arbitrary  constant.  Thus  in 
equation  (5)  K  may  only  be  reduced  to  a  nonsingular  matrix  by  specifying  val¬ 
ues  for  u(t ! )  and/or  u(t2)  so  that  the  variations  of  one  or  both  of  these 
quantities  vanish.  The  essence  of  the  discussion  which  follows  is  not  changed 
if,  in  the  sequel,  it  is  assumed  that  u(t^)  has  been  specified.  This  is  known 
as  a  ’geometric'  or  'imposed'  constraint.  Because  5Uj  =  SuCtj)  =  0  multiplies 
the  first  row  of  K  in  equation  (5),  this  row  is  effectively  removed  from  the 
formulation.  Since  the  remaining  variations  are  arbitrary  the  final  set  of 
equations  to  be  solved  is  then: 
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i  =  2,3. ..n 


(9) 


l  Kijuj  “  Fi  “  Kilul  * 
j=2 


where  =  u(t^)  is  the  specified  value  and  n  is  the  dimension  of  K.  Whether 
these  equations  derive  from  a  global  power  series  approximation  or  from  one 
based  on  finite  elements,  one  may  readily  verify  that  as  n  is  increased  their 
solutions  do  indeed  converge  to  the  exact  solution  of  the  corresponding  two 
point  time-boundary  value  problem.  Should  one  wish  a  solution  to  an  initial 
value  problem,  however,  equation  (4)  must  be  used  instead  of  equation  (5).  In 
this  case,  specifying  values  for  u(t^)  and  u(t^)  cause  5U]_  and  6U2  to  vanish 
thereby  deleting  the  first  two  equations  of  this  set.  The  resulting  system  of 
equations  to  be  solved  is  thus: 
n 

l  (Kij-Bij)Uj  =  Fi  -  (Kil-Bil)U1-(Ki2-Bi2)U2  >  1  =  3»4 . n  0°) 

j=3 

In  all  cases  attempted  to  date,  solutions  to  equations  (10)  have  been  observed 
to  converge  to  the  exact  solution  if  these  equations  are  derived  using  a 
global  power  series  approximation  but  not  if  they  are  formulated  by  finite 
elements.  An  example  of  this  anomaly  will  be  given  in  the  next  section.  As 
the  only  difference  between  equations  (4)  and  (5)  is  a  subtraction  of  B  in  the 
former,  and  in  as  much  as  convergence  is  achieved  when  equation  (4)  derives 
from  a  power  series  approximation,  one  suspects  that  it  is  the  finite  element 
representation  of  the  matrix  B  which  is  somehow  at  fault.  It  is  therefore  of 
interest  to  know  in  more  detail  just  how  the  subtraction  of  B  is  supposed  to 
affect  the  coefficient  matrix  of  the  system. 

In  contrast  to  the  matrix  K,  the  matrix  K-B  must  tend  to  be  singular  of 
degeneracy  two  -  no  constraints  having  been  assumed  a  priori.  Thus  when  u(t^) 
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is  specified  and  the  first  row  of  K-B  is  deleted,  the  remaining  equations 
still  must  possess  one  degeneracy  in  the  limit  as  the  number  of  basis  func¬ 
tions  becomes  infinite.  Thus  the  effect  of  subtracting  B  must  be  to  free  the 
natural  boundary  condition  at  t2  (inherent  in  equation  (5))  and  to  introduce. a 
degeneracy.  This  remaining  degeneracy  can  only  be  removed  by  specifying  the 

value  of  u(t)  at  a  time  other  than  or  a  value  for  u,  resulting  in  the  dele¬ 
tion  of  another  row  of  K-B. 

ANOMALOUS  BEHAVIOR  OF  FINITE  ELEMENT  FORMULATIONS 

The  degree  to  which  the  subtraction  of  the  matrix  B  from  K  can  both  free 
the  natural  boundary  condition  at  t2  and  introduce  a  degeneracy  differs  with 
the  type  of  approximation  employed.  When  global  power  series  approximations 
are  usee  the  B  matrix  is  quite  full  and  the  subtraction  affects  many  rows  of 
K.  When  locally  defined  Hermite  polynomials  are  used,  however,  B  is  very 
sparse  and  in  fact  contains  only  two  non-zero  components.  Moreover,  one  of 
these  appears  in  the  first  row  of  B  which  is  deleted  when  u(tL)  is  specified. 

In  this  case  freeing  the  natural  boundary  condition  and  introducing  a  degener¬ 
acy  depends  on  the  subtraction  from  a  single  component  of  K.  Even  though  both 
effects  may  actually  be  produced  in  the  limit  as  the  number  of  elements 
becomes  infinite,  the  degree  to  which  they  are  approximated  for  any  finite 
number  of  elements  is  evidently  insufficient  and  the  solutions  do  not  converge 
to  the  correct  result.  This  is  exemplified  in  Figure  1.  The  problem  repres¬ 
ented  is  that  of  a  free  oscillator  of  unit  mass  and  stiffness,  subject  to  the 
prescribed  initial  constraints  of  zero  displacement  and  unit  velocity.  For 
this  case,  equation  (2a)  reads: 
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(11) 


/  (u6u-u6u)dt  -  u6u|  -  0 
0  0 

or  simply, 

J  (u-f-u)6udt  =  0  *  (12) 

0 


The  finite  element  results  of  Figure  1  were  obtained  using  piecewise  cubic 
Hermite  polynomials.  (Higher  ordered  Hermite  polynomials  yield  similar 
results.)  It  Is  observed  that  the  solutions  tend  to  diminish  from  the  exact 
solution,  sin(t),  as  the  number  of  elements  is  increased.  Using  only  two  fin 
ite  elements  the  finite  element  matrix  formulation  (equation  (4))  for  this 
problem  is  as  follows: 

0  =  <5Ut[K-B]  U  =  [5Ui  <SU2  6U3  6U4  6U5  6U6]  • 


Jkll 

ki  2 

k13 

k14 

0 

1 

0  i 

I 

0  -1 

0 

0 

0 

1 

ol 

1 

1  Ui  1 

Jk21 

k22 

k23 

k24 

0 

0  1 

1 

0 

0 

0 

0 

0 

ol 

1 

1  02 

1  1 

jk31 

k32 

k33+kU 

k34+k12 

k13 

k14  1 

0 

0 

0 

0 

0 

ol 

1  U3  1 

1  k41 

k42 

k43+k21 

k44+k22 

k23 

k24  1 

0 

0 

0 

0 

0 

ol 

1 

1 

1  U4  1 

1  1 

|0 

0 

k31 

k32 

k33 

k34  1 

I 

0 

0 

0 

0 

0 

1 

1 1 

1 

1  u5 

io 

0 

k4l 

k42 

k43 

1 

k44 1 

1 

0 

0 

0 

0 

0 

ol 

1 

1  o6 

1  1 

(13) 


*Note  that  Eq .  (12)  would  also  result  from  application  of  the  Galerkin 
procedure,  implying  that  the  Galerkin  method  has  physical  justification  for 
problems  In  dynamics. 
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Using  expression  (8),  the  element  matrix  k  is  calculated  in  terms  of  the 
element  length  h  as: 


6  13h  1  llh2  9h  6  13h2  1 


5h  35  10  210  70  5h  420  10 


2h  h3  13h2  1  h3  h 


15  105  420  10  140  30 


6  13h  llh2  1 


SYMM 


5h  35  210  10 


2h  h3 


15  105 


Since  is  specified  the  first  row  of  K  -  B  is  deleted.  As  the  subtraction 
of  B  only  affects  one  row  of  the  reduced  system,  the  only  way  in  which  a 
degeneracy  can  be  introduced  is  for  the  next  to  last  row  to  join  the  space 
defined  by  the  rows  remaining.  Thus  rows  two  through  six  in  equation  (13) 
ideally  would  become  linearly  dependent.  This  dependency  among  rows  must  be 
quite  general  as  specification  of  any  other  of  the  Uj_  must  remove  it. 

One  suspects  that  a  simple  subtraction  of  unity  from  K55  in  equation  (13) 
may  not  do  the  best  job  of  introducing  a  degeneracy  or  of  freeing  the  natural 
boundary  condition  at  t2  =  71  •  One  can  gain  some  idea  of  how  ’close1  this 
subtraction  brings  the  fifth  row  into  the  space  of  rows  2,3,4  and  6  by 
comparing  it  with  its  projection  onto  this  space.  Substituting  tt/2  for  h,  the 
fifth  row  of  equation  (13)  calculates  to  be: 


*A11  mathematics  herein  were  performed  using  the  MACSYMA  (Project  MAC’s 
SYmbolic  MAnipulat ion)  system  developed  by  the  Mathlab  Group  of  the  MIT 
Laboratory  for  Computer  Science. 
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[0.0  0.0  -0.96590326  -0.17637194  0.180505097  -0.970755175] 

whereas  its  projection  is: 

[7.8587183E-3  -8.5978979E-3  -0.974496335 

-0.184380835  0.172642875  -0.9617834]  . 

Further  calculations  show  that  if  the  interval  of  solution  remains  fixed  and 
the  number  of  finite  elements  is  allowed  to  increase,  closer  agreement  between 
the  next  to  last  row  vector  and  its  projection  is  observed  but  this  is  not 
accompanied  by  a  convergence  of  the  solution  vector  toward  the  exact  solution 
to  the  problem.  While  the  exact  reasons  for  this  instability  are  not  known  it 
is  apparent  that  the  rate  at  which  the  next  to  last  row  tends  to  become  depen¬ 
dent  is  important.  It  stands  to  reason,  therefore,  that  should  one  invoke  the 
limit  condition  without  actually  proceeding  to  the  limit,  a  convergent 
sequence  may  result  and  Indeed  this  proves  to  be  the  case. 

Asserting  that  the  row  vectors  two  through  six  are  linearly  dependent 
allows  the  fifth  row  (equation)  of  equations  (13)  to  be  replaced  by  a  linear 
combination  of  the  others.  For  example,  let 

R5  =  a2R2  +  a3R3  +  a4R4  +  a6R6  (14) 

where  denotes  the  i^h  row  of  R  —  B.  After  imposing  the  second  initial 
constraint,  U2  =  1 ,  equations  (13)  can  be  written: 

6U3R3  •  U  +  6U4R4  •  U  +  6U5(a2R2+a3R3+a4R4+a5R5)  •  U  +  5U5R5  •  U  =  0  (15) 

Since  all  variations  in  equation  (15)  are  arbitrary,  there  results  the 
following  system  of  equations  for  solution: 

0  =  R3  •  U  =  R2  •  U  =  R4  •  U  =  R6  •  U  (16) 
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Thus  the  second  equation  (row)  which  was  originally  deleted  through  the  speci¬ 
fication  of  U2,  is  brought  back,  into  the  formulation  in  place  of  the  fifth  in 
a  logical  and  consistent  manner.  Equations  (16)  are  the  same  set  as  would 
result  from  following  the  procedure  of  Argyris  and  Scharpf.  These  authors, 
however,  started  with  Hamilton’s  principle  which  requires  that  =  6U5  =  0. 

This  would  delete  the  first  and  fifth  equations  from  the  set.  Further  speci¬ 
fication  of  U2  should  then  delete  the  second  equation  as  well,  overspecifying 
the  problem.  Argyris  and  Scharpf ^  allow  this  equation  to  remain  without  jus¬ 
tification.  Moreover,  no  explanation  is  given  as  to  why  6U5  should  vanish  as 
U5  is  never  specified  in  an  initial  value  problem.  All  of  these  inconsisten¬ 
cies  derive  from  the  fact  that  Hamilton’s  principle  corresponds  only  to  bound¬ 
ary  value  problems  -  never  to  initial  value  problems. 

In  summary,  the  work  of  this  section  shows  that  Hamilton’s  law  of  varying 
action,  unlike  Hamilton’s  principle,  is  an  unconstrained  variational  statement 
permitting  the  introduction  of  arbitrary  constraints  including  data  ordinarily 
given  for  initial  value  problems.  When  piecewise  Hermite  cubic  polynomials 
are  used  as  a  basis  for  a  finite  element  formulation,  the  singular  state  of 
the  resulting  coefficient  matrix  in  the  limit  justifies  retention  of  the 
second  equation  of  the  system  in  preference  to  the  next  to  last  when  typical 
initial  values  for  displacement  and  velocity  are  specified.  Following  this 
procedure,  convergent  solutions  are  then  obtained  for  the  problem  of  the  free 
oscillator  considered  in  this  section.  These  results  are  presented  in  Table  I 
for  formulations  based  on  one,  two,  and  six  finite  elements. 

7 Argyris,  J.  H. ,  and  Scharpf,  D.  W. ,  "Finite  Elements  in  Time  and  Space," 
Nuclear  Engineering  and  Design,  10,  1969,  456-464. 
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TABLE  I.  SOLUTIONS  TO  FREE  OSCILLATOR  PROBLEM  (DISPLACEMENT/VELOCITY) 


0  <  t  <  TT 


6t/n 

One  Element 

Two  Elements 

Six  Elements 

Exact 

Solution 

0 

0.0* 

0.0* 

0.0* 

0.0 

1.0* 

1.0* 

1.0* 

1.0 

1 

0.49978005 

0.5 

0.86602547 

0.86602541 

2 

0.86564452 

0.86602541 

0.50000025 

0.5 

3 

0.97817298 

0.99956036 

1.0 

2.02985945E-4 

4.4572957E-7 

0.0 

4 

0.86564496 

0.86602541 

-0.49999948 

-0.5 

5 

0.499780823 

0.5 

-0.86602502 

0.86602541 

.6 

0.0166090783 

3.9845105E-4 

8.9120273E-7 

0.0 

-1.00079414 

-1.00000946 

-0.99999999 

-1.0 

*  Imposed  values. 
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APPLICATIONS 


Example  1.  Linear  Oscillator  Subjected  to  Discontinuous  Forces 
A  linear  oscillator  of  unit  mass  and  stiffness  is  subjected  to  a  force 
f(t).  Two  cases  are  considered: 

(a)  f(t)  =  H(t-l/2) 

(b)  f(t)  =  6 (t-1/2) 

H  and  6  are  the  Heaviside  and  Dirac  functions  respectively  and  for  either  of 
these  cases  equation  (2)  reads: 

1 2  •  •  •  1 2 

/  {u6u  +  ( f (t )-u)6u} dt  -  u5u  |  =  0  (17) 

For  case  (a)  four  finite  elements  of  equal  length  are  used  to  approximate  u(t) 
over  the  solution  interval  (0,2) •  The  element  polynomial  shape  function  is 
Hermite  cubic  and  an  element  length  of  one  half  takes  advantage  of  the  specif¬ 
ic  shape  of  the  forcing  function.  Table  II  compares  the  calculated  displace¬ 
ments  and  velocities  with  those  computed  from  the  exact  solution. 

In  case  (b)  a  discontinuity  in  velocity  can  be  expected  in  the  solution. 
As  the  use  of  cubic  shape  functions  enforces  continuity  of  velocity  through¬ 
out,  a  better  solution  might  be  expected  when  linear  shape  functions  are 
employed.  Table  III  compares  the  exact  solution  on  the  interval  (0,1)  with 
that  obtained  using  ten  such  elements  of  equal  length. 

The  two  problems  considered  in  this  example  demonstrate  the  manner  in 
which  the  type  of  element  and  its  points  of  attachment  (i.e.,  the  ’nodes’  or 
’grid  points’)  may  be  varied  to  suit  specified  transient  events. 
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TABLE  II.  SOLUTION  TO  u  +  u  =  H(t-l/2) 


0  <  t  <  2.0 


t 

Computed 

Displacement  Velocity 

Displacement 

Exact 

Velocity 

0.0 

0.0* 

1.0* 

0.0 

1.0 

0.5 

0.47932149 

0.87708716 

0.47942555 

0.877582565 

1.0 

0.96370936 

1.0199163 

0.96388844 

1.01972786  ■ 

1.5 

1.45700388 

0.91238744 

1.45719267 

0.91220819 

2.0 

1.83836447 

0.5805616 

1.83856024 

0.58134814 

*Iraposed  values. 

TABLE  III.  SOLUTION  TO  u  +  u  =  6(t-l/2) 


0  <  t  <  1 

t 

Computed  Displacement 

Exact  Displacement 

0.0 

0.0* 

0.0 

0.1 

0.1* 

0.099833416 

0.2 

0.199001664 

0.19866933 

0.3 

0.296016622 

0.295520213 

0.4 

0.390076343 

0.38941834 

0.5 

0.58007539 

0.57925896 

0.6 

0.76428335 

0.76331182 

0.7 

0.94086118 

0.93973791 

0.8 

1.10804607 

1.10677443 

0.9 

1.26416892 

1.26275246 

1 .0 

1.40767112 

1.40611348 

* Imposed  values. 
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Example  2.  Response  of  a  Beam  to  a  Moving  Mass 

A  concentrated  mass  is  assumed  to  move  at  constant  velocity  v  along  the 
length  of  a  uniform  Euler  beam,  simply  supported  at  each  of  its  ends  and 
having  zero  displacement  and  velocity  at  t  =  0.  Under  suitable  definitions 
for  k  and  m,  the  representative  equations  may  be  written:^ 

yiv  +  ky  +  f(x,t)  =  0 

y(0,t)  =  y" (0  ,t )  =  y(l,t)  =  y" ( 1 ,t)  =  y(x,0)  =  y(x,0)  =  0  (18) 

The  function  f(x,t)  consists  of  a  sum  of  inertial  terms: 

•  — 

f(x,t)  =  m(y  +  2vy’  +  g  +  v2y") 5(x-vt )  (19) 

where  g  denotes  the  gravitational  constant  and  6  is  the  Dirac  function.  This 
problem  is  particularly  interesting  in  that  the  conventional  use  of  piecewise 
cubic  shape  functions  to  discretize  the  space  variable  only,  introduces  forces 
which  are  discontinuous  functions  of  time  into  the  resulting  ordinary  differ¬ 
ential  equations.  These  discontinuities  are  associated  with  the  beam  curva¬ 
ture  load  term  appearing  in  the  expression  (19).  Since  the  piecewise  cubic 
polynomials  are  discontinuous  in  the  second  derivative  at  the  element  attach¬ 
ments,  the  term  mv2y"<5  (x-vt )  -  when  multiplied  by  the  shape  function  a(x)  and 
integrated  over  the  element  length  -  will  produce  functions  of  time  which  are 
discontinuous  whenever  the  moving  mass  arrives  at  any  point  of  attachment. 
Clearly  these  discontinuities  have  nothing  to  do  with  the  physics  of  the  prob¬ 
lem  anc  are  certain  to  invite  trouble  when  one  attempts  to  numerically  inte- 

l^Simkins,  T.  E. ,  "Unconstrained  Variational  Statements  for  Initial  and 
Boundary-Value  Problems,"  AIAA  Journal,  Vol.  16,  No.  6,  June  1978,  pp. 
559-563. 
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grate  the  time  dependent  equations  via  established  algorithms.  It  is  possi¬ 
ble,  of  course,  to  use  shape  functions  of  higher  degree  to  discretize  the 
space  variable  thus  eliminating  the  discontinuities  at  the  onset  but  this  is 
hardly  consistent  with  the  finite  element  method  which  should  permit  the  use 
of  even  linear  shape  functions  if  need  be.  One  is  tempted  to  somehow 
’smooth1  these  discontinuities,  yet  this  should  not  be  done  in  a  purely 
arbitrary  fashion.  Integrating  the  effects  of  these  forces  throughout  the 
time  domain  through  the  use  of  Hamilton’s  law  of  varying  action  provides  a 
consistent  way  to  handle  this  problem. 

While  it  is  possible  to  handle  the  space  and  time  finite  element 
discretizations  in  one  operation,  the  amount  of  computation  and  computer 
programming  tend  to  become  inordinately  large.  Moreover,  there  exist  any 
number  of  finite  element  codes  (e.g.  NASTRAN)  which  can  quickly  accomplish 
much  of  the  space  discretization.  It  seems  more  efficient,  therefore,  to 
apply  the  finite  element  method  in  two  steps,  by  first  discretizing  the  space 
variable  and  then  applying  Hamilton’s  law  to  the  resulting  system  of  ordinary 
differential  equations  in  time.  For  the  case  at  hand,  the  differential 
equations  governing  the  motion  of  the  i*-^1  beam  element  turn  out  to  be: 

A  •  •  A  A  A 

(p  +  mc^)u  +  mc2u  +  (q  +  mc3)u  +  mga(vt)=  0  (20) 

p  and  q  are  proportional  to  the  usual  mass  and  stiffness  matrices  for  beam 
elements  and  have  been  evaluated  many  times  in  the  literature.  Here  all  of 
the  beam  elements  are  of  the  same  length  £,  and  the  displacement  within  the 
itVl  element  is  interpolated  from  u*(t),  a  vector  of  end  point  displacements 
and  velocities,  i.e., 
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y(x,t)  =  a^(5^)u^(t) 

0  <  S1  <  1 


(21) 


where  £i(x)  =  x/£  -  ( i— 1 ) ,  a  nondiraenslonal  element  coordinate. 

The  c  matrices  In  equations  (20)  correspond  to  transverse,  Coriolis,  and- 
centrifugal  accelerations  respectively  and  are  defined  for  the  i1-'1  element  as 
follows: 

ex  =  a(^)  aTCC^Ix-vt 

C2  =  2va(5^)a,T(5^) lx=vt  (22) 

C3  =  v2a(5i)a"T(5i) lx=vt 

It  is  noted  that  C3  will  be  discontinuous  at  =  0  and  ■  1.  The  function 
m  takes  on  the  value  of  m  only  when  the  concentrated  mass  lies  within  the  ith 

A 

element,  otherwise  m  is  zero. 

The  element  equations  (20)  are  combined  in  the  usual  way  to  form  N  equa¬ 
tions  of  motion  for  the  combined  structure.  Symbolically: 

**  • 

M(t)U  +  C(t)U  +  K(t)U  =  F(t )  (23) 

Each  of  the  matrices  in  equation  (23)  can  be  viewed  as  a  conventional  matrix 
of  constant  coefficients  plus  a  time  variant  set  of  components  which  are 
active  in  a  band  along  its  main  diagonal  as  the  moving  mass  traverses  the  beam 
in  time.  For  this  system  of  equations  Hamilton’s  law  of  varying  action  can  be 
written : 

l  I  {/t2{6U1MijUJ  +  6Ui[(Mij-Cij)Uj-KijUj  +  F^ldt  -  611^1^  |2  }  =  0  (24) 

i=l  j=l  t!  tl 
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It  is  interesting  to  observe  the  accuracy  of  solution  which  can  be 
obtained  from  equation  (24)  using  only  two  finite  elements  in  space  and  two  in 
time.  A  formulation  using  two  elements  in  space  results  in  a  system  of  N=4 
ordinary  differential  equations  in  time  once  the  geometric  support  constraints 
have  been  applied.  A  two  element  formulation  of  these  four  equations  for  the 
time  domain,  followed  by  the  application  of  all  initial  constraints  in  the 
manner  summarized  in  Section  5,  gives  a  final  system  of  sixteen  linear 
algebraic  equations  for  solution.  Figure  2  compares  this  solution  with  the 
experimental  results  of  Ayre ,  Jacobsen,  and  Hsu1^  and  a  conventional  finite 
element  solution  using  three  elements  in  the  space  domain  followed  by  a 
time-integration  of  the  equations  (28)  by  Hamming’s  predictor-corrector 

algorithm. The  mass  velocity  in  this  case  is  v  =  v*/2,  - v*  being  the 

lowest  velocity  to  cause  resonance  when  the  load  is  a  moving  weight  only  and 
the  magnitude  assigned  to  the  moving  mass  is  25%  of  the  total  mass  of  the 
beam.  (Other  parametric  values  are  the  same  as  those  in  reference  16.)  The 
displacements  have  been  normalized  with  respect  to  the  maximum  deflection 
produced  if  the  weight  was  applied  statically  at  midspan  and  L  is  the  total 
beam  length.  In  particular  one  notes  that  the  conventional  solution  obtained 
via  three  finite  elements  in  space  only,  produces  non-physical  discontinuities 
in  the  slope  of  the  solution  curve  at  vt/L  =  1/3,  2/3.  (The  continuous  data 

l^Ayre,  R.  S. ,  Jacobsen,  L.  S.,  and  Hsu,  C.  S. ,  "Transverse  Vibration  of  One 
and  of  Two  Space  Beams  Under  the  Action  of  a  Moving  Mass  Load,"  Proceedings 
of  First  National  Congress  on  Applied  Mechanics,  June  1951. 

^Ralston  and  Wilf,  Mathematical  Methods  for  Digital  Computers,  Wiley  and 
Sons,  NY,  London,  1960,  pp.  95-109. 
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for  generating  this  curve  is  obtained  by  interpolating  the  solution  to 
equation  (23)  using  equation  (21).)  No  discontinuities  of  this  sort  can  arise 
when  finite  elements  in  space  and  time  are  employed.  Improved  agreement  with 
the  experimental  results  is  also  observed. 
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Figure  1.  Divergent  Finite  Element  Solutions  to  Free  Oscillator  Problems 
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Figure  2,  Displacement  of  Beam  at  Location  of  Moving  Mass. 
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